<!DOCTYPE html>
<html>
<head><meta charset="utf-8" />
<title>Maxwell3ND1femrate</title><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>

<style type="text/css">
/* Overrides of notebook CSS for static HTML export */
body {
  overflow: visible;
  padding: 8px;
}
div#notebook {
  overflow: visible;
  border-top: none;
}@media print {
  div.cell {
    display: block;
    page-break-inside: avoid;
  } 
  div.output_wrapper { 
    display: block;
    page-break-inside: avoid; 
  }
  div.output { 
    display: block;
    page-break-inside: avoid; 
  }
}
</style>

<!-- Custom stylesheet, it must be in the parent directory as the html file -->
<link rel="stylesheet" type="text/css" media="all" href="../doc.css" />
<link rel="stylesheet" type="text/css" media="all" href="doc.css" />

<!-- Loading mathjax macro -->
<!-- Load mathjax -->
    <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS_HTML"></script>
    <!-- MathJax configuration -->
    <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
        tex2jax: {
            inlineMath: [ ['$','$'], ["\\(","\\)"] ],
            displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
            processEscapes: true,
            processEnvironments: true
        },
        // Center justify equations in code and markdown cells. Elsewhere
        // we use CSS to left justify single line equations in code cells.
        displayAlign: 'center',
        "HTML-CSS": {
            styles: {'.MathJax_Display': {"margin": 0}},
            linebreaks: { automatic: true }
        }
    });
    </script>
    <!-- End of mathjax configuration --></head>
<body>
  <div tabindex="-1" id="notebook" class="border-box-sizing">
    <div class="container" id="notebook-container">

<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="Linear-Edge-Element-for-Maxwell-Equations-in-3D">Linear Edge Element for Maxwell Equations in 3D<a class="anchor-link" href="#Linear-Edge-Element-for-Maxwell-Equations-in-3D">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>This example is to show the linear edge element approximation of the electric field of the time harmonic Maxwell equation.</p>
\begin{align}
\nabla \times (\mu^{-1}\nabla \times  u) - \omega^2 \varepsilon \, u &amp;= J  \quad  \text{ in } \quad \Omega,  \\
                                  n \times u &amp;= n \times g_D  \quad  \text{ on } \quad \Gamma_D,\\
                    n \times (\mu^{-1}\nabla \times  u) &amp;= n \times g_N  \quad  \text{ on } \quad \Gamma_N.
\end{align}<p>based on the weak formulation</p>
$$(\mu^{-1}\nabla \times  u, \nabla \times  v) - (\omega^2\varepsilon u,v) = (J,v) - \langle n \times g_N,v \rangle_{\Gamma_N}.$$
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p><strong>Reference</strong></p>
<ul>
<li><a href="http://www.math.uci.edu/~chenlong/226/FEMMaxwell.pdf">Finite Element Methods for Maxwell Equations</a></li>
<li><a href="http://www.math.uci.edu/~chenlong/226/codeMaxwell.pdf">Programming of Finite Element Methods for Maxwell Equations</a></li>
</ul>
<p><strong>Subroutines</strong>:</p>

<pre><code>- Maxwell1
- cubeMaxwell1
- femMaxwell3
- Maxwell1femrate

</code></pre>
<p>The method is implemented in <code>Maxwell1</code> subroutine and tested in <code>cubeMaxwell1</code>. Together with other elements (ND0,ND1,ND2), <code>femMaxwell3</code> provides a concise interface to solve Maxwell equation. The ND1 element is tested in <code>Maxwell1femrate</code>. This doc is based on <code>Maxwell1femrate</code>.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Data-Structure">Data Structure<a class="anchor-link" href="#Data-Structure">&#182;</a></h2><p>Use the function</p>

<pre><code>[elem2dof,edge,elem2edgeSign] = dof3edge(elem);

</code></pre>
<p>to construct the pointer from element index to edge index. Read</p>
<p>&lt;dof3edgedoc.html Dof on Edges in Three Dimensions&gt; for details.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[2]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">node</span> <span class="p">=</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">];</span>
<span class="n">elem</span> <span class="p">=</span> <span class="p">[</span><span class="mi">1</span> <span class="mi">2</span> <span class="mi">3</span> <span class="mi">4</span><span class="p">];</span>
<span class="n">localEdge</span> <span class="p">=</span> <span class="p">[</span><span class="mi">1</span> <span class="mi">2</span><span class="p">;</span> <span class="mi">1</span> <span class="mi">3</span><span class="p">;</span> <span class="mi">1</span> <span class="mi">4</span><span class="p">;</span> <span class="mi">2</span> <span class="mi">3</span><span class="p">;</span> <span class="mi">2</span> <span class="mi">4</span><span class="p">;</span> <span class="mi">3</span> <span class="mi">4</span><span class="p">];</span>
<span class="n">imatlab_export_fig</span><span class="p">(</span><span class="s">&#39;print-png&#39;</span><span class="p">)</span>  <span class="c">% Static png figures.</span>
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Units&#39;</span><span class="p">,</span><span class="s">&#39;normal&#39;</span><span class="p">);</span> 
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Position&#39;</span><span class="p">,[</span><span class="mf">0.25</span><span class="p">,</span><span class="mf">0.25</span><span class="p">,</span><span class="mf">0.25</span><span class="p">,</span><span class="mf">0.25</span><span class="p">]);</span>
<span class="n">showmesh3</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">view</span><span class="p">(</span><span class="o">-</span><span class="mi">72</span><span class="p">,</span><span class="mi">9</span><span class="p">);</span>
<span class="n">findnode3</span><span class="p">(</span><span class="n">node</span><span class="p">);</span>
<span class="n">findedge</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">localEdge</span><span class="p">,</span><span class="s">&#39;all&#39;</span><span class="p">,</span><span class="s">&#39;vec&#39;</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>The six dofs associated to edges in a tetrahedron is sorted in the ordering <code>[1 2; 1 3; 1 4; 2 3; 2 4; 3 4]</code>. Here <code>[1 2 3 4]</code> are local indices of vertices.</p>
<p>Globally we use ascend ordering for each element and thus the orientation of the edge is consistent. No need of <code>elem2edgeSign</code>. Read <a href="../mesh/sc3doc.html">Simplicial complex in three dimensions</a> for more discussion of indexing, ordering and orientation.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Local-Bases">Local Bases<a class="anchor-link" href="#Local-Bases">&#182;</a></h2><p>Suppose <code>[i,j]</code> is the kth edge and <code>i&lt;j</code>. The basis is given by</p>
$$ \phi _k = \lambda_i\nabla \lambda_j - \lambda_j \nabla \lambda_i,\qquad
   \nabla \times \phi_k = 2\nabla \lambda_i \times \nabla \lambda_j.$$<p>Inside one tetrahedron, the 6 bases functions along with their curl
corresponding to 6 local edges <code>[1 2; 1 3; 1 4; 2 3; 2 4; 3 4]</code> are</p>
$$ \phi_1 = \lambda_1\nabla\lambda_2 - \lambda_2\nabla\lambda_1,\qquad
   \nabla \times \phi_1 = 2\nabla\lambda_1\times \nabla\lambda_2,$$$$ \phi_2 = \lambda_1\nabla\lambda_3 - \lambda_3\nabla\lambda_1,\qquad
   \nabla \times \phi_2 = 2\nabla\lambda_1\times \nabla\lambda_3,$$$$ \phi_3 = \lambda_1\nabla\lambda_4 - \lambda_4\nabla\lambda_1,\qquad
   \nabla \times \phi_3 = 2\nabla\lambda_1\times \nabla\lambda_4,$$$$ \phi_4 = \lambda_2\nabla\lambda_3 - \lambda_3\nabla\lambda_2,\qquad
   \nabla \times \phi_4 = 2\nabla\lambda_2\times \nabla\lambda_3,$$$$ \phi_5 = \lambda_2\nabla\lambda_4 - \lambda_4\nabla\lambda_2,\qquad
   \nabla \times \phi_5 = 2\nabla\lambda_2\times \nabla\lambda_4,$$$$ \phi_6 = \lambda_3\nabla\lambda_4 - \lambda_4\nabla\lambda_3,\qquad
   \nabla \times \phi_6 = 2\nabla\lambda_3\times \nabla\lambda_4.$$<p>The additional 6 bases for the second family are:</p>
$$ \psi_k = \lambda_i\nabla \lambda_j + \lambda_j \nabla \lambda_i,\qquad
   \nabla \times \psi_k = 0.$$$$ \psi_1 = \lambda_1\nabla\lambda_2 + \lambda_2\nabla\lambda_1,$$$$ \psi_2 = \lambda_1\nabla\lambda_3 + \lambda_3\nabla\lambda_1,$$$$ \psi_3 = \lambda_1\nabla\lambda_4 + \lambda_4\nabla\lambda_1,$$$$ \psi_4 = \lambda_2\nabla\lambda_3 + \lambda_3\nabla\lambda_2,$$$$ \psi_5 = \lambda_2\nabla\lambda_4 + \lambda_4\nabla\lambda_2,$$$$ \psi_6 = \lambda_3\nabla\lambda_4 + \lambda_4\nabla\lambda_3.$$
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Degree-of-freedoms">Degree of freedoms<a class="anchor-link" href="#Degree-of-freedoms">&#182;</a></h2><p>Suppose <code>[i,j]</code> is the kth edge and <code>i&lt;j</code>. The corresponding degree of freedom is</p>
$$l_k (v) = \int_{e_k} v\cdot t \, {\rm d}s \approx \frac{1}{2}[v(i)+v(j)]\cdot e_{k}.$$<p>It is dual to the basis $\{\phi_k\}$ in the sense that</p>
$$l_{\ell}(\phi _k) = \delta_{k,\ell}.$$<p>The additional 6 degree of freedoms are:</p>
$$l_k^1 (v) = 3\int_{e_k} v\cdot t(\lambda _i - \lambda_j) \, {\rm d}s  \approx \frac{1}{2}[v(i) - v(j)]\cdot e_{k}.$$
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Dirichlet-boundary-condition">Dirichlet boundary condition<a class="anchor-link" href="#Dirichlet-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[3]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">%% Setting</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">cubemesh</span><span class="p">([</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span><span class="mi">1</span><span class="p">);</span>
<span class="n">mesh</span> <span class="p">=</span> <span class="n">struct</span><span class="p">(</span><span class="s">&#39;node&#39;</span><span class="p">,</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;elem&#39;</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">option</span><span class="p">.</span><span class="n">L0</span> <span class="p">=</span> <span class="mi">1</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">maxIt</span> <span class="p">=</span> <span class="mi">4</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">elemType</span> <span class="p">=</span> <span class="s">&#39;ND1&#39;</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">printlevel</span> <span class="p">=</span> <span class="mi">1</span><span class="p">;</span>

<span class="c">%% Dirichlet boundary condition.</span>
<span class="n">fprintf</span><span class="p">(</span><span class="s">&#39;Dirichlet boundary conditions. \n&#39;</span><span class="p">);</span>    
<span class="n">pde</span> <span class="p">=</span> <span class="n">Maxwelldata2</span><span class="p">;</span>
<span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary3</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Dirichlet&#39;</span><span class="p">);</span>
<span class="n">femMaxwell3</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Dirichlet boundary conditions. 
#dof:     1208, Direct solver  0.1 
#dof:     8368, Direct solver 0.23 
Conjugate Gradient Method using HX preconditioner 
#dof:    62048,   #nnz:  1372268,   iter: 37,   err = 6.5415e-09,   time =  1.6 s
Conjugate Gradient Method using HX preconditioner 
#dof:   477376,   #nnz: 11791756,   iter: 36,   err = 9.7991e-09,   time =   11 s
Table: Error
 #Dof        h        ||u-u_h||    ||Du-Du_h||   ||DuI-Du_h|| ||uI-u_h||_{max}

  1208   2.500e-01   5.52165e-02   3.96067e-01   7.95505e-02   3.10119e-03
  8368   1.250e-01   1.36258e-02   1.97200e-01   4.02371e-02   4.33068e-04
 62048   6.250e-02   3.39303e-03   9.84696e-02   2.02922e-02   5.78792e-05
477376   3.125e-02   8.47174e-04   4.92119e-02   1.01987e-02   7.53487e-06

Table: CPU time
 #Dof    Assemble     Solve      Error      Mesh    

  1208   9.00e-02   1.00e-01   4.00e-02   1.00e-02
  8368   1.20e-01   2.30e-01   7.00e-02   1.00e-02
 62048   7.90e-01   1.56e+00   2.20e-01   5.00e-02
477376   6.15e+00   1.07e+01   1.71e+00   0.00e+00

</pre>
</div>
</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Pure-Neumann-boundary-condition">Pure Neumann boundary condition<a class="anchor-link" href="#Pure-Neumann-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[5]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">%% Pure Neumann boundary condition.</span>
<span class="n">fprintf</span><span class="p">(</span><span class="s">&#39;Neumann boundary condition. \n&#39;</span><span class="p">);</span>
<span class="n">option</span><span class="p">.</span><span class="n">plotflag</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">Maxwelldata2</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary3</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Neumann&#39;</span><span class="p">);</span>
<span class="n">femMaxwell3</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Neumann boundary condition. 
#dof:     1208, Direct solver 0.06 
#dof:     8368, Direct solver 0.35 
Conjugate Gradient Method using HX preconditioner 
#dof:    62048,   #nnz:  1628576,   iter: 43,   err = 8.7403e-09,   time =    2 s
Conjugate Gradient Method using HX preconditioner 
#dof:   477376,   #nnz: 12838720,   iter: 43,   err = 9.4520e-09,   time =   13 s
Table: Error
 #Dof        h        ||u-u_h||    ||Du-Du_h||   ||DuI-Du_h|| ||uI-u_h||_{max}

  1208   2.500e-01   2.41877e-02   3.90295e-01   9.93550e-02   1.88901e-02
  8368   1.250e-01   6.28974e-03   1.96199e-01   4.40961e-02   2.94393e-03
 62048   6.250e-02   1.59327e-03   9.82971e-02   2.10068e-02   4.34979e-04
477376   3.125e-02   4.00195e-04   4.91812e-02   1.03323e-02   6.24133e-05

Table: CPU time
 #Dof    Assemble     Solve      Error      Mesh    

  1208   6.00e-02   6.00e-02   1.00e-02   0.00e+00
  8368   1.80e-01   3.50e-01   4.00e-02   1.00e-02
 62048   7.00e-01   1.95e+00   2.20e-01   4.00e-02
477376   5.47e+00   1.35e+01   1.69e+00   0.00e+00

</pre>
</div>
</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Conclusion">Conclusion<a class="anchor-link" href="#Conclusion">&#182;</a></h2><p>The H(curl)-norm is still 1st order but the L2-norm is improved to 2nd order.</p>
<p>MGCG using HX preconditioner converges uniformly in all cases.</p>

</div>
</div>
</div>
    </div>
  </div>
</body>

 


</html>
